import numpy as np
from scipy.ndimage import convolve
import os
import glob
from skimage import io, img_as_float
from skimage.transform import resize
import matplotlib.pyplot as plt
from time import time

import cv2
import time
from skimage import io
from skimage.filters import gaussian
from skimage.morphology import disk, dilation

# 计时开始
start_time = time.time()

#%% 加载图像函数
def LoadingImg(folder, prefix, start_idx, end_idx, ext, step):
    images = []
    for i in range(start_idx, end_idx + 1, step):
        filename = f"{folder}{prefix}{i}{ext}"
        img = io.imread(filename)  # 读取图像
        img_gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)  # 转换为灰度图
        images.append(img_gray)
    H, W = images[0].shape
    nbframe = len(images)
    return images, H, W, nbframe

#%% 生成对焦地图函数
def In_focus_Maps(images, sigma, threshold_for_quantile, se_size):
    infocus_maps = []
    for img in images:
        # 计算高斯梯度
        grad_x = cv2.Sobel(gaussian(img, sigma), cv2.CV_64F, 1, 0, ksize=3)
        grad_y = cv2.Sobel(gaussian(img, sigma), cv2.CV_64F, 0, 1, ksize=3)
        # 计算梯度的幅度
        grad_magnitude = np.sqrt(grad_x**2 + grad_y**2)
        
        # 使用分位数设定阈值
        threshold = np.quantile(grad_magnitude, threshold_for_quantile)
        threshold_map = grad_magnitude > threshold
        
        # 膨胀处理
        se = disk(se_size)
        infocus_map = dilation(threshold_map, se)
        infocus_maps.append(infocus_map)
    return np.array(infocus_maps)

#%% 生成金字塔函数
def StackGP(images, a):
    G = []
    for img in images:
        gp = img
        pyramid = [gp]
        while gp.shape[0] > 1 and gp.shape[1] > 1:
            gp = cv2.pyrDown(gp.astype(np.uint8), dstsize=(int(gp.shape[1] * a), int(gp.shape[0] * a)))
            pyramid.append(gp)
        G.append(pyramid)
    return {"G": G}

def StackLP(images, a):
    L = []
    for img in images:
        gp = img
        pyramid = []
        while gp.shape[0] > 1 and gp.shape[1] > 1:
            next_gp = cv2.pyrDown(gp, dstsize=(int(gp.shape[1] * a), int(gp.shape[0] * a)))
            lp = cv2.subtract(gp, cv2.pyrUp(next_gp, dstsize=(gp.shape[1], gp.shape[0])))
            pyramid.append(lp)
            gp = next_gp
        pyramid.append(gp)  # 最后一层高斯金字塔
        L.append(pyramid)
    return {"L": L}

#%% 融合金字塔
def PyramidFusion(IBWLP, a):
    fused_pyramid = []
    num_levels = len(IBWLP['L'][0])
    for level in range(num_levels):
        fused_layer = np.zeros_like(IBWLP['L'][0][level])
        for i in range(len(IBWLP['L'])):
            fused_layer += IBWLP['L'][i][level]
        fused_layer = fused_layer / len(IBWLP['L'])
        fused_pyramid.append(fused_layer)
    return {"G": fused_pyramid}

#%% 金字塔逆变换
def ILP(FusedLP, a):
    gp = FusedLP['G'][-1]
    for i in range(len(FusedLP['G']) - 2, -1, -1):
        gp = cv2.pyrUp(gp, dstsize=(FusedLP['G'][i].shape[1], FusedLP['G'][i].shape[0]))
        gp = cv2.add(gp, FusedLP['G'][i])
    return {"G": [gp]}

#%% 主程序
# 加载图像
folder = 'Demo_from_Helicon_Focus/'  # 图像文件夹
prefix = 'DSC0'  # 图像前缀
start_idx = 1  # 起始索引
end_idx = 3  # 结束索引
ext = '.jpg'  # 图像格式
step = 1  # 步长
IBW, H, W, nbframe = LoadingImg(folder, prefix, start_idx, end_idx, ext, step)

# 生成对焦地图
sigma = 0.5  # 高斯梯度参数
threshold_for_quantile = 0.8  # 分位数阈值
se_size = 10  # 膨胀核大小
InfocusMaps = In_focus_Maps(IBW, sigma, threshold_for_quantile, se_size)

# 生成灰度金字塔
a = 0.5
InfocusGP = StackGP(InfocusMaps, a)
IBWLP = StackLP(IBW, a)

# 预选对焦像素
for i in range(len(InfocusGP['G'][0])):
    for j in range(len(IBWLP['L'])):
        IBWLP['L'][j][i] = IBWLP['L'][j][i] * InfocusGP['G'][j][i]

# 融合金字塔
FusedLP = PyramidFusion(IBWLP, a)

# 生成全对焦灰度图像
FusedLP = ILP(FusedLP, a)
Iref = FusedLP['G'][0][:H, :W]

# 显示融合图像
cv2.imshow('Fused Image', Iref.astype(np.uint8))
cv2.waitKey(0)
cv2.destroyAllWindows()

# 计时结束
end_time = time.time()
print(f"运行时间: {end_time - start_time:.2f} 秒")




# Expand.m
# def kernel(a):
#     """
#     Generates a 1D Gaussian kernel based on the parameter 'a'.
#     The full kernel (3D) is created later.
#     """
#     size = int(4 * a + 1)  # Define kernel size
#     x = np.linspace(-size // 2, size // 2, size)
#     gaussian_1d = np.exp(-(x**2) / (2 * a**2))
#     gaussian_1d /= gaussian_1d.sum()  # Normalize
#     return np.outer(gaussian_1d, gaussian_1d)[:, :, np.newaxis]  # Convert to 3D kernel


# # def expand(gl1, a):
# #     """
# #     Expand function that reconstructs a Gaussian pyramid layer.

# #     Parameters:
# #     - gl1: Input layer (numpy array) to be expanded
# #     - a: Parameter for the Gaussian kernel

# #     Returns:
# #     - gl0: Expanded output layer
# #     """
# #     # Generate the Gaussian kernel
# #     w = kernel(a)

# #     # Get the size of the current layer
# #     rl1, cl1, nbframe = gl1.shape

# #     # Initialize an empty array for the expanded layer
# #     gl0 = np.zeros(((rl1 - 1) * 2 + 1, (cl1 - 1) * 2 + 1, nbframe))

# #     # Place the input pixels into the expanded grid
# #     gl0[::2, ::2, :] = gl1

# #     # Convolve with the Gaussian kernel (3D convolution across frames)
# #     for i in range(nbframe):
# #         gl0[:, :, i] = 4 * convolve(gl0[:, :, i], w[:, :, 0], mode='constant', cval=0.0)

# #     return gl0


# # GoGFilter.m
# def gog_filter(image, sigma):
#     """
#     Computes the magnitude and argument of the Gaussian gradient of the image.

#     Parameters:
#     - image: Input 2D numpy array (grayscale image).
#     - sigma: Gaussian smoothing parameter (typically less than 1).

#     Returns:
#     - gog_mag: Magnitude of the Gaussian gradient.
#     - gog_arg: Argument (angle) of the Gaussian gradient.
#     """
#     gog_x, gog_y = gog_filter_xy(image, sigma)
#     gog_mag = np.sqrt(gog_x ** 2 + gog_y ** 2)
#     gog_arg = np.arctan2(gog_y, gog_x)
#     return gog_mag, gog_arg


# def gog_filter_xy(image, sigma):
#     """
#     Computes the Gaussian gradient of the image in X and Y directions.

#     Parameters:
#     - image: Input 2D numpy array (grayscale image).
#     - sigma: Gaussian smoothing parameter.

#     Returns:
#     - gog_x: Gradient along the X direction.
#     - gog_y: Gradient along the Y direction.
#     """
#     # Define kernel size
#     x = np.arange(-2, 3)  # Equivalent to -2:2 in MATLAB
#     y = x
#     X, Y = np.meshgrid(x, y)

#     # Gaussian gradient kernels
#     gog_x_kernel = np.exp(-(X ** 2 + Y ** 2) / (2 * sigma ** 2)) * (-X) / (2 * np.pi * sigma ** 4)
#     gog_y_kernel = np.exp(-(X ** 2 + Y ** 2) / (2 * sigma ** 2)) * (-Y) / (2 * np.pi * sigma ** 4)

#     # Convolve the input image with the Gaussian gradient kernels
#     gog_x = convolve(image, gog_x_kernel, mode='constant', cval=0.0)
#     gog_y = convolve(image, gog_y_kernel, mode='constant', cval=0.0)

#     return gog_x, gog_y


# # ILP.m
# def ilp(pyr, a):
#     """
#     逆构建拉普拉斯金字塔。

#     参数:
#     - pyr: 包含高斯金字塔和拉普拉斯金字塔的字典，结构如下：
#            pyr['G'] - 高斯金字塔列表
#            pyr['L'] - 拉普拉斯金字塔列表
#     - a: 高斯核参数

#     返回:
#     - pyr: 更新后的金字塔，完成逆构建。
#     """
#     lvl_max = len(pyr['G'])  # 高斯金字塔的最大层数

#     # 从倒数第二层开始，逐层合并拉普拉斯金字塔
#     for i in range(lvl_max - 2, -1, -1):  
#         print(f"Collapse pyramid at level : {i}")
#         pyr['G'][i] = expand(pyr['G'][i + 1], a) + pyr['L'][i]
#         pyr['G'][i + 1] = None  # 清空下一层，释放空间
    
#     return pyr


# def expand(gl1, a):
#     """
#     扩展高斯金字塔层到上一层的尺寸。

#     参数:
#     - gl1: 输入的高斯金字塔层（numpy数组）。
#     - a: 高斯核参数。

#     返回:
#     - gl0: 扩展后的图像层。
#     """
#     from scipy.ndimage import convolve

#     # 生成高斯核
#     w = kernel(a)

#     # 获取输入图像尺寸
#     rl1, cl1 = gl1.shape[:2]
#     gl0 = np.zeros(((rl1 - 1) * 2 + 1, (cl1 - 1) * 2 + 1))

#     # 将输入像素放置到扩展图像的偶数位置
#     gl0[::2, ::2] = gl1

#     # 使用高斯核进行卷积，放大图像并平滑
#     gl0 = 4 * convolve(gl0, w, mode='constant', cval=0.0)
#     return gl0



# # In _focus_Maps.m
# def in_focus_maps(IBW, sigma, threshold_for_quantile, se_size):
#     """
#     计算图像中聚焦像素的地图。

#     参数:
#     - IBW: 输入的灰度图像序列（numpy数组，形状为 HxWxnbframe）。
#     - sigma: 高斯滤波的平滑参数。
#     - threshold_for_quantile: 分位数阈值，用于二值化操作。
#     - se_size: 结构元素的大小，用于膨胀操作。

#     返回:
#     - InfocusMaps: 膨胀后的聚焦像素地图。
#     """
#     H, W, nbframe = IBW.shape  # 获取图像的尺寸

#     ## 生成边缘地图
#     print("生成边缘地图...")
#     edges_maps, _ = gog_filter(IBW, sigma)  # 高斯梯度的幅值
#     weighting_edges_maps = edges_maps**2 / np.sum(edges_maps**2, axis=2, keepdims=True)

#     ## 基于分位数生成二值掩码
#     print("生成二值掩码...")
#     binary_masks = np.zeros((H, W, nbframe), dtype=bool)
#     for i in range(nbframe):
#         print(f"生成二值掩码 {i+1}/{nbframe}")
#         flattened = weighting_edges_maps[:, :, i].reshape(-1)
#         threshold = mquantiles(flattened, threshold_for_quantile)[0]  # 计算分位数阈值
#         binary_masks[:, :, i] = weighting_edges_maps[:, :, i] > threshold

#     ## 扩展聚焦像素地图
#     print("扩展聚焦像素地图...")
#     se = diamond(se_size)  # 生成菱形结构元素
#     infocus_maps = np.zeros_like(binary_masks)
#     for i in range(nbframe):
#         infocus_maps[:, :, i] = binary_dilation(binary_masks[:, :, i], se)

#     return infocus_maps


# def gog_filter(IBW, sigma):
#     """
#     计算图像的高斯梯度幅值和方向。

#     参数:
#     - IBW: 输入的灰度图像序列。
#     - sigma: 高斯平滑参数。

#     返回:
#     - gog_mag: 高斯梯度的幅值。
#     - gog_arg: 高斯梯度的方向。
#     """
#     H, W, nbframe = IBW.shape
#     gog_mag = np.zeros((H, W, nbframe))
#     gog_arg = np.zeros((H, W, nbframe))
#     for i in range(nbframe):
#         edges_mag, edges_arg = gog_filter_single(IBW[:, :, i], sigma)
#         gog_mag[:, :, i] = edges_mag
#         gog_arg[:, :, i] = edges_arg
#     return gog_mag, gog_arg


# def gog_filter_single(image, sigma):
#     """
#     计算单个图像的高斯梯度幅值和方向。

#     参数:
#     - image: 输入的2D图像。
#     - sigma: 高斯平滑参数。

#     返回:
#     - gog_mag: 高斯梯度的幅值。
#     - gog_arg: 高斯梯度的方向。
#     """
#     # 生成高斯核的梯度
#     x = np.arange(-2, 3)
#     y = x
#     X, Y = np.meshgrid(x, y)
#     gog_x = np.exp(-(X**2 + Y**2) / (2 * sigma**2)) * (-X) / (2 * np.pi * sigma**4)
#     gog_y = np.exp(-(X**2 + Y**2) / (2 * sigma**2)) * (-Y) / (2 * np.pi * sigma**4)

#     # 卷积计算梯度
#     Ix = convolve(image, gog_x, mode='constant', cval=0.0)
#     Iy = convolve(image, gog_y, mode='constant', cval=0.0)

#     # 计算幅值和方向
#     gog_mag = np.sqrt(Ix**2 + Iy**2)
#     gog_arg = np.arctan2(Iy, Ix)

#     return gog_mag, gog_arg



# # Kernel.m
# def Kernel(a):
#     """
#     该函数计算一个用于扩展和缩减金字塔图像的核。
    
#     参数:
#     a (float): 用于核计算的参数。
    
#     返回:
#     numpy.ndarray: 计算出的二维核。
#     """
#     # 计算一维核
#     w1d = np.array([1/4 - a/2, 1/4, a, 1/4, 1/4 - a/2])
    
#     # 计算二维核
#     w2d = np.outer(w1d, w1d)
    
#     return w2d


# # Loadinglmg.m


# def LoadingImg(directory, imgname, nostart, noend, imgext, imgsize):
#     """
#     该函数用于加载指定目录下从nostart到noend的图像，并根据需要调整图像大小。

#     参数:
#     directory (str): 图像所在的目录路径。
#     imgname (str): 图像文件的前缀名称。
#     nostart (int): 起始图像编号。
#     noend (int): 结束图像编号。
#     imgext (str): 图像的扩展名（例如'.jpg'）。
#     imgsize (float): 图像缩放因子。如果为0或1，则不调整图像大小。

#     返回:
#     IRGB (ndarray): 存储RGB图像的数组。
#     IBW (ndarray): 存储灰度图像的数组。
#     H (int): 图像的高度。
#     W (int): 图像的宽度。
#     nbframe (int): 加载的图像帧数。
#     """
#     nbframe = noend - nostart + 1
#     # 获取图像的尺寸信息
#     imgpath = os.path.join(directory, f"{imgname}{nostart}{imgext}")
#     imginfo = cv2.imread(imgpath)
#     H, W, _ = imginfo.shape  # 获取原始图像的尺寸
    
#     # 根据imgsize调整图像大小
#     if imgsize != 0 and imgsize != 1:
#         H = int(np.ceil(H * imgsize))
#         W = int(np.ceil(W * imgsize))
#         ResizeImg = True
#     else:
#         ResizeImg = False
    
#     # 初始化存储空间
#     IRGB = np.zeros((H, W, 3 * nbframe), dtype=np.uint8)
#     IBW = np.zeros((H, W, nbframe), dtype=np.float32)
    
#     idx = 0
#     # 加载图像并计算灰度图像
#     for i in range(nbframe):
#         print(f'Loading Img n° {i+1}/{nbframe}')
#         filename = os.path.join(directory, f"{imgname}{nostart + i}{imgext}")
        
#         # 读取并调整图像大小（如果需要）
#         img = cv2.imread(filename)
#         if ResizeImg:
#             img = cv2.resize(img, (W, H))
        
#         IRGB[:, :, idx:idx+3] = img
#         IBW[:, :, i] = RGB2BW(img)
#         idx += 3
    
#     return IRGB, IBW, H, W, nbframe


# def RGB2BW(rgb):
#     """
#     该函数将RGB图像转换为灰度图像。

#     参数:
#     rgb (ndarray): 输入的RGB图像。

#     返回:
#     gray (ndarray): 输出的灰度图像。
#     """
#     if rgb.ndim != 3 or rgb.shape[2] != 3:
#         raise ValueError("输入的图像不是RGB图像")
    
#     # 将RGB图像转换为灰度图像
#     gray = np.sum(rgb, axis=2) / 3.0
#     return gray



# # PaddZeros.m
# def PaddZeros(I1):
#     """
#     该函数将图像I1用零填充，使得图像的大小变为离其最近的2的幂次方加1。

#     参数:
#     I1 (ndarray): 输入图像，通常是一个三维数组（高度 x 宽度 x 帧数）。

#     返回:
#     I2 (ndarray): 填充后的图像。
#     LvlMax (int): 图像的最大级别（2的幂次方的指数）。
#     """
#     row, col, Nbframe = I1.shape
    
#     # 计算行和列的最近2的幂次方的指数
#     Nrow = np.ceil(np.log2(row)).astype(int)
#     Ncol = np.ceil(np.log2(col)).astype(int)
    
#     # 选择最大的指数值
#     N = max(Nrow, Ncol)
    
#     # 计算最终尺寸
#     S = (2 ** N) + 1
    
#     # 创建一个全零的输出数组
#     I2 = np.zeros((S, S, Nbframe), dtype=I1.dtype)
    
#     # 将原始图像复制到输出图像的左上角
#     I2[:row, :col, :] = I1
    
#     # 返回填充后的图像和最大级别
#     LvlMax = N
#     return I2, LvlMax


# # Pyramid.m
# def Pyramid(Level, type):
#     """
#     根据金字塔的类型和高度（级别）定义金字塔的数据。

#     参数:
#     Level (int): 金字塔的层数。
#     type (str): 金字塔的类型，可以是 'Laplacian' 或 'Gaussian'。

#     返回:
#     PyrDef (dict): 包含金字塔数据的字典。
#     """
#     if type == 'Laplacian':
#         PyrDef = {'G': [None] * Level, 'L': [None] * Level}
#     elif type == 'Gaussian':
#         PyrDef = {'G': [None] * Level}
#     else:
#         raise ValueError(f"The type {type} is unknown")
    
#     return PyrDef


# # PyramidFusion.m
# import numpy as np
# import cv2
# from scipy.ndimage import convolve,  standard_deviation
# from scipy.stats import entropy

# def PyramidFusion(Pyrs, a):
#     """
#     该函数根据局部区域能量计算金字塔的融合，每一层都使用局部能量标准，顶层使用熵和标准差标准。

#     参数:
#     Pyrs (dict): 包含金字塔数据的字典，必须包含字段 'G' 和 'L'。
#     a (float): 用于生成核的参数。

#     返回:
#     FusPyr (dict): 融合后的金字塔。
#     """
#     # 获取金字塔的层数
#     LvlMax = len(Pyrs['G'])
    
#     # 初始化融合金字塔
#     FusPyr = Pyramid(LvlMax, 'Laplacian')
    
#     # 应用局部区域能量标准（除顶层外的所有层）
#     w = Kernel(a)
#     for lvl in range(LvlMax - 1):
#         print(f'Fusing pyramids at level {lvl + 1}')
        
#         # 计算局部区域能量
#         LocalRegionEnergy = np.abs(convolve(Pyrs['L'][lvl] ** 2, w, mode='nearest'))
        
#         # 找到最大能量对应的索引
#         idx = np.argmax(LocalRegionEnergy, axis=2)
        
#         # 获取图像的行列数
#         r, c, _ = Pyrs['L'][lvl].shape
#         N = r * c
        
#         # 根据索引重新排列金字塔的层
#         idx = [i + (idx.flatten() - 1) * N for i in range(N)]
#         FusPyr['L'][lvl] = Pyrs['L'][lvl].reshape(-1)[idx].reshape(r, c)

#     # 应用最大熵和标准差标准（对于顶层）
#     idxE = np.argmax(entropy(Pyrs['G'][LvlMax - 1]), axis=2)
#     idxD = np.argmax(standard_deviation(Pyrs['G'][LvlMax - 1]), axis=2)
    
#     # 获取图像的行列数
#     r, c, _ = Pyrs['G'][LvlMax - 1].shape
#     N = r * c
    
#     # 重新排列顶层图像
#     idx = [i + (idxE.flatten() - 1) * N for i in range(N)]
#     GendE = Pyrs['G'][LvlMax - 1].reshape(-1)[idx].reshape(r, c)
    
#     idx = [i + (idxD.flatten() - 1) * N for i in range(N)]
#     GendD = Pyrs['G'][LvlMax - 1].reshape(-1)[idx].reshape(r, c)
    
#     # 计算融合后的顶层图像
#     FusPyr['G'][LvlMax - 1] = (GendE + GendD) / 2
    
#     return FusPyr


# def Pyramid(Level, type):
#     """
#     根据金字塔的类型和高度（级别）定义金字塔的数据。

#     参数:
#     Level (int): 金字塔的层数。
#     type (str): 金字塔的类型，可以是 'Laplacian' 或 'Gaussian'。

#     返回:
#     PyrDef (dict): 包含金字塔数据的字典。
#     """
#     if type == 'Laplacian':
#         PyrDef = {'G': [None] * Level, 'L': [None] * Level}
#     elif type == 'Gaussian':
#         PyrDef = {'G': [None] * Level}
#     else:
#         raise ValueError(f"The type {type} is unknown")
    
#     return PyrDef

# # Reduce.m

# def Reduce(Gl0, a):
#     """
#     将图像 Gl0 通过高斯模糊和下采样操作缩小一半。

#     参数:
#     Gl0 (ndarray): 输入的三维图像数组 (高度 x 宽度 x 帧数)。
#     a (float): 参数 a，用于模糊操作时生成高斯核。

#     返回:
#     Gl1 (ndarray): 下采样后的图像数组。
#     """
#     # 生成高斯模糊核
#     w = Kernel(a)
    
#     # 使用卷积对图像进行模糊处理
#     Gl0_blurred = convolve(Gl0, w, mode='nearest')
    
#     # 下采样：每隔2个像素选取一个
#     Gl1 = Gl0_blurred[::2, ::2, :]
    
#     return Gl1



# # StackGP.m
# def StackGP(I, a):
#     """
#     计算图像 I 的高斯金字塔。

#     参数:
#     I (ndarray): 输入图像数组 (高度 x 宽度 x 帧数)。
#     a (float): 高斯模糊核的参数，用于生成高斯-like 核。

#     返回:
#     GP (dict): 包含高斯金字塔的字典，包含多个金字塔层次。
#     """
#     # 填充图像到最近的2的幂次方加1
#     I, Levels = PaddZeros(I)
    
#     # 初始化高斯金字塔
#     GP = Pyramid(Levels, 'Gaussian')
    
#     # 生成金字塔每一层
#     for k in range(Levels):
#         print(f'Gaussian Pyramid level {k+1}/{Levels}')
        
#         if k == 0:
#             GP['G'][0] = I
#         else:
#             GP['G'][k] = Reduce(GP['G'][k-1], a)
    
#     return GP

# def PaddZeros(I1):
#     """
#     该函数将图像 I1 用零填充，使得图像的大小变为离其最近的2的幂次方加1。

#     参数:
#     I1 (ndarray): 输入图像。

#     返回:
#     I2 (ndarray): 填充后的图像。
#     LvlMax (int): 图像的最大级别（2的幂次方的指数）。
#     """
#     row, col, Nbframe = I1.shape
#     Nrow = np.ceil(np.log2(row)).astype(int)
#     Ncol = np.ceil(np.log2(col)).astype(int)
#     N = max(Nrow, Ncol)
#     S = (2 ** N) + 1
#     I2 = np.zeros((S, S, Nbframe), dtype=I1.dtype)
#     I2[:row, :col, :] = I1
#     return I2, N

# def Pyramid(Level, type):
#     """
#     根据金字塔的类型和高度（级别）定义金字塔的数据。

#     参数:
#     Level (int): 金字塔的层数。
#     type (str): 金字塔的类型，可以是 'Laplacian' 或 'Gaussian'。

#     返回:
#     PyrDef (dict): 包含金字塔数据的字典。
#     """
#     if type == 'Laplacian':
#         PyrDef = {'G': [None] * Level, 'L': [None] * Level}
#     elif type == 'Gaussian':
#         PyrDef = {'G': [None] * Level}
#     else:
#         raise ValueError(f"The type {type} is unknown")
    
#     return PyrDef

# def Reduce(Gl0, a):
#     """
#     将图像 Gl0 通过高斯模糊和下采样操作缩小一半。

#     参数:
#     Gl0 (ndarray): 输入的三维图像数组 (高度 x 宽度 x 帧数)。
#     a (float): 参数 a，用于模糊操作时生成高斯核。

#     返回:
#     Gl1 (ndarray): 下采样后的图像数组。
#     """
#     # 生成高斯模糊核
#     w = Kernel(a)
    
#     # 使用卷积对图像进行模糊处理
#     Gl0_blurred = convolve(Gl0, w, mode='nearest')
    
#     # 下采样：每隔2个像素选取一个
#     Gl1 = Gl0_blurred[::2, ::2, :]
    
#     return Gl1


# # StackLP.m


# def StackLP(I, a):
#     """
#     计算图像 I 的拉普拉斯金字塔。

#     参数:
#     I (ndarray): 输入图像数组 (高度 x 宽度 x 帧数)。
#     a (float): 用于生成高斯模糊核的参数。

#     返回:
#     LP (dict): 包含拉普拉斯金字塔的字典，包含多个金字塔层次。
#     """
#     # 填充图像到最近的2的幂次方加1
#     I, Levels = PaddZeros(I)
    
#     # 初始化拉普拉斯金字塔
#     LP = Pyramid(Levels, 'Laplacian')
    
#     # 生成金字塔每一层
#     for k in range(Levels):
#         print(f'Pyramids level {k+1}/{Levels}')
        
#         if k == 0:
#             LP['G'][0] = I
#         else:
#             LP['G'][k-1] = None
        
#         # 生成下一层
#         LP['G'][k+1] = Reduce(LP['G'][k], a)
        
#         # 计算拉普拉斯层
#         LP['L'][k] = LP['G'][k] - Expand(LP['G'][k+1], a)
    
#     LP['G'][Levels] = None
    
#     return LP

# def Expand(Gl0, a):
#     """
#     将图像 Gl0 放大一倍，使用高斯核进行模糊。

#     参数:
#     Gl0 (ndarray): 输入的三维图像数组。
#     a (float): 参数 a，用于模糊操作时生成高斯核。

#     返回:
#     Gl1 (ndarray): 放大后的图像。
#     """
#     # 生成高斯模糊核
#     w = Kernel(a)
    
#     # 使用卷积对图像进行模糊处理
#     Gl0_blurred = convolve(Gl0, w, mode='nearest')
    
#     # 放大图像：插值放大
#     Gl1 = np.repeat(Gl0_blurred, 2, axis=0)
#     Gl1 = np.repeat(Gl1, 2, axis=1)
    
#     return Gl1

# def PaddZeros(I1):
#     """
#     该函数将图像 I1 用零填充，使得图像的大小变为离其最近的2的幂次方加1。

#     参数:
#     I1 (ndarray): 输入图像。

#     返回:
#     I2 (ndarray): 填充后的图像。
#     LvlMax (int): 图像的最大级别（2的幂次方的指数）。
#     """
#     row, col, Nbframe = I1.shape
#     Nrow = np.ceil(np.log2(row)).astype(int)
#     Ncol = np.ceil(np.log2(col)).astype(int)
#     N = max(Nrow, Ncol)
#     S = (2 ** N) + 1
#     I2 = np.zeros((S, S, Nbframe), dtype=I1.dtype)
#     I2[:row, :col, :] = I1
#     return I2, N

# def Pyramid(Level, type):
#     """
#     根据金字塔的类型和高度（级别）定义金字塔的数据。

#     参数:
#     Level (int): 金字塔的层数。
#     type (str): 金字塔的类型，可以是 'Laplacian' 或 'Gaussian'。

#     返回:
#     PyrDef (dict): 包含金字塔数据的字典。
#     """
#     if type == 'Laplacian':
#         PyrDef = {'G': [None] * Level, 'L': [None] * Level}
#     elif type == 'Gaussian':
#         PyrDef = {'G': [None] * Level}
#     else:
#         raise ValueError(f"The type {type} is unknown")
    
#     return PyrDef

# def Reduce(Gl0, a):
#     """
#     将图像 Gl0 通过高斯模糊和下采样操作缩小一半。

#     参数:
#     Gl0 (ndarray): 输入的三维图像数组 (高度 x 宽度 x 帧数)。
#     a (float): 参数 a，用于模糊操作时生成高斯核。

#     返回:
#     Gl1 (ndarray): 下采样后的图像数组。
#     """
#     # 生成高斯模糊核
#     w = Kernel(a)
    
#     # 使用卷积对图像进行模糊处理
#     Gl0_blurred = convolve(Gl0, w, mode='nearest')
    
#     # 下采样：每隔2个像素选取一个
#     Gl1 = Gl0_blurred[::2, ::2, :]
    
#     return Gl1



